/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */



/*! \file FinitieDifferenceDelta.h
    \brief Multidimensional generalization of finite difference Laplacian operators.

	Multidimensional generalization of the classical Laplacian difference operator
	\f[
			\delta_{x\overline{x}} u_i=\frac{u_{i-1}+2u_i-u_{i-1}}{h^2}
	\f]
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/
#if !defined(_FD_Delta_)
#define _FD_Delta_
#include <src/miscellaneous/OS_basic.h>
#include <src/math/Operators/FiniteDifferenceOperators.h>
#include <src/math/BoundaryCondition/ClassicalGraphDefinition.h>
#include <src/algorithm/Evolution_Algorithm/CFL_traits.h>
#include <src/algorithm/Evolution_Algorithm/Predictor_traits.h>
#include <src/data/OS_Multi_Array.h>
#pragma warning (disable: 4355)
/*! 
    \brief Multidimensional Laplacian operator.

	Multidimensional Laplacian operators, i.e.
	\f[
			\Delta_{x\overline{x}} u_{i_1,..,i_D}=\sum_{j=1,..,D}\frac{u_{i_1,..,i_j+1,..,i_D}-2u_{i_1,..,i_D}+u_{i_1,..,i_j-1,..,i_D}}{h_i^2}
	\f]
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/

template <class data, std::size_t D, class V = smart_ptr<data>>
class Delta 
	: public FiniteDifferenceOperator<data,data,V>
 {
public :
		OS_STATIC_CHECK(Dim > 1);
		enum:std::size_t												{Order = 2};
		typedef OTSCentered												Centered_Type;
		typedef typename data::Boundary_type							Boundary_type;
		
		typedef FiniteDifferenceOperator								Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef Delta													type;
		typedef typename data::value_type								value_type;

//! return " Delta u".
		static std::string get_name(){return " Delta u ";};
		Delta( ) : operator_interface_(*this,make_3<base_value_type>(0.,0.,2.)) {};
		Delta(const Delta & data) : operator_interface_(*this),  Base(data) {};
		Delta(const weight_set_type & data) :  operator_interface_(*this,make_3<base_value_type>(0.,0.,2.)), Base(data) {};

		Delta(const Base_Cast &first) : operator_interface_(*this), Base(first)
		{};
		virtual ~Delta(){};
		Delta & operator = (const Delta & from) {
			Base::operator =(from);
			return *this;
		};
		virtual void update() {
			operator_interface_.update();
			Delta_1_.set_state_variable( get_state_variable() ) ;
		};
		//! Optimized code for finite difference Laplacian operators.
		virtual result_type operator()(typename argument_type from) const 
		{
			typedef argument_base_type::const_iterator const_iterator;
			typedef result_type::iterator iterator;
			result_type result;
			size_type size = operator_interface_.size();
			if (size >0) {
				result = ( (from>>1) - 2.*from + (from<<1))*(size*size);
				const_iterator beg = from.begin();
				const_iterator end = from.end();
				iterator beg_r = result.begin();
				for (;beg != end; ++beg, ++beg_r) {
					*beg_r += Delta_1_(*beg);
				};
			};
			return result;
		};
		virtual void build_graph() {
			OS_STATIC_CHECK(Dim == 2); // to implement later when needed.
			typedef weight_set_type::const_iterator1 const_iterator1;
			typedef weight_set_type::const_iterator2 const_iterator2;

			int DIM (state_variable_type::value_type::Dim)  ;

			indice_type size1 = get_state_variable()->size(1);
			indice_type size2 = get_state_variable()->size(2);
			graph_.clear();
			for (indice_type i = 0; i < size1;++i)
				for (indice_type j = 0; j < size2;++j)
				{
					indice_type indice_i = i+j*size1;
					indice_type indice = i+(std::min(j+1, size2-1) )*size1; //up
					graph_.insert(std::make_pair(indice_i,indice));
					indice = i+(std::max(j-1, 0) )*size1; //down
					graph_.insert(std::make_pair(indice_i,indice));
					indice = std::max(i-1, 0) +j*size1; //left
					graph_.insert(std::make_pair(indice_i,indice));
					indice = std::min(i+1, size1-1) +j*size1; //left
					graph_.insert(std::make_pair(indice_i,indice));
				}
		};

		//! Factory method.
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};
//		typedef typename smart_ptr<typename data::reference> data_state_variable_type;
		operator_interface<Delta,Order,Centered_Type,Boundary_type> operator_interface_;
		Delta<typename data::value_type, Dim-1, V> Delta_1_;
};

/*! 
    \brief Laplacian on vectors
	If the input has protected access out of bounds (as OS_Vector), the Laplacian is not defined.
	We suppose in this case that the Laplacian is defined over each component let \f$\{u^k\}_{k=1,..,N}\f$ be a vector of object, then
	this specialization computes
	\f[
			\Delta_{x\overline{x}} u^k_{i_1,..,i_D}=\sum_{j=1,..,D}\frac{u^k_{i_1,..,i_j+1,..,i_D}-2u^k_{i_1,..,i_D}+u^k_{i_1,..,i_j-1,..,i_D}}{h_i^2}
	\f]
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/
template <class data, std::size_t D, class V >
class Delta <OS_Vector<data>, D, V>
	: public FiniteDifferenceOperator<OS_Vector<data>,OS_Vector<data>,V>
 {
public :
		OS_STATIC_CHECK(Dim == 3);
		enum:std::size_t												{Order = 2};
		typedef OTSCentered												Centered_Type;
		typedef typename data::Boundary_type							Boundary_type;
		
		typedef FiniteDifferenceOperator								Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef typename Base::graph_type                               graph_type;
		typedef Delta													type;
		typedef typename data::value_type								value_type;

//! return " Delta u".
		static std::string get_name(){return " {Delta u}^k ";};
		type( ) {};
		type(const type & data) : Base(data) {};
		type(const Base_Cast &first) : operator_interface_(*this), Base(first)
		{};
		virtual ~Delta(){};
		type & operator = (const type & from) {
			Base::operator =(from);
			return *this;
		};
		virtual void update() {
			Delta_1_.set_state_variable( get_state_variable() );
		};
		//! Optimized code for finite difference Laplacian operators.
		virtual result_type operator()(typename argument_type from) const 
		{
			typedef argument_base_type::const_iterator const_iterator;
			typedef result_type::iterator iterator;
			if (from.capacity() == 0)
				return result_type();
			result_type result(from.size_list());
			const_iterator beg = from.begin();
			const_iterator end = from.end();
			iterator beg_r = result.begin();
			for (;beg != end; ++beg, ++beg_r) {
				OS_Array_Copy(*beg_r, Delta_1_(*beg));
			};
			return result;
		};
		virtual void build_graph() {
			Delta_1_.build_graph();
		};
		virtual const graph_type& get_graph() const {
			return Delta_1_.get_graph();
		};

		//! Factory method.
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,type>(this);
		};
		Delta<data, Dim-1,V> Delta_1_;
};


/*! \brief One dimensional Laplacian operator specialization.

	One dimensional Laplacian operators, i.e.
	\f[
			\delta_{x\overline{x}} u_{i}=\frac{u_{i-1}-2u_{i}+u_{i+1}}{h^2}
	\f]
	The geometry is the unit cube \f$[0,1]^N\f$, where N is given by data::Dim
*/

template <class data, class V>
class Delta<data,1,V>  
	: public FiniteDifferenceOperator<data,data,V>
{
public :
		OS_STATIC_CHECK(data::Dim == 1);
		enum:std::size_t												{Order = 2};
		typedef OTSCentered												Centered_Type;
		typedef typename data::Boundary_type							Boundary_type;
		
		typedef FiniteDifferenceOperator								Base;
		typedef  typename Base::Base_Cast								Base_Cast;
		typedef Delta													type;

//! return " Delta u".
		static std::string get_name(){return " D_xx u ";};
		Delta( ) : operator_interface_(*this,make_3<base_value_type>(0.,0.,2.)) {};
		Delta(const Delta & data) : operator_interface_(*this),  Base(data) {};
		Delta(const weight_set_type & data) :  operator_interface_(*this,make_3<base_value_type>(0.,0.,2.)), Base(data) {};

		Delta(const Base_Cast &first) : operator_interface_(*this), Base(first)
		{};
		virtual ~Delta(){};
		Delta & operator = (const Delta & from) {
			Base::operator =(from);
			return *this;
		};
		virtual void update() {
			operator_interface_.update();
		};
		//! Optimized code for finite difference Laplacian operators.
		virtual result_type operator()(typename argument_type from) const 
		{
			size_type size = operator_interface_.size();
			if (size >0)
				return (argument_type)( (from>>1) - 2.*from + (from<<1))*(size*size);
			return result_type();
		};
		//! Factory method.
		virtual smart_ptr_type Create_Instance() const {
			return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
		};
		operator_interface<Delta,Order,Centered_Type,Boundary_type> operator_interface_;
};

/*! \brief Bind operator of Laplacian difference operator Delta_traits

	This is the time Bind operator of Laplacian difference operator Delta_traits, that is Delta_t(t,u)=Delta(u)
*/
template <class data, std::size_t Dim, typename Time = OS_double >
struct Delta_t : public BindTimeOperator<Delta<data, Dim> ,Time > {

	typedef Delta_t														time_operator_type;
	typedef BindTimeOperator<time_operator_type ,Time >					Base;
	Delta_t(){};
};



#endif